Imaging apparatus and method

ABSTRACT

A method of imaging a sample, the method comprising: (a) irradiating a sample with a pulse of electromagnetic radiation, said pulse having a plurality of frequencies in the range from 25 GHz to 100 THz; (b) determining a first parameter related to the amplitude of the radiation, which is either reflected from and/or transmitted by the sample, in the time domain; (c) calculating the value of the first parameter at a first time value relative to the value of the first parameter at a second time value which coincides with a physical feature of the dataset of the first parameter with respect to time; and (d) generating an image by plotting the value calculated in step (c) for different points of the sample.

BACKGROUND OF THE INVENTION

The present invention relates to the field of imaging samples withradiation in the infra-red (IR) and Terahertz frequency range andspecifically using radiation in the higher Gigahertz (GHz) and theTerahertz (THz) frequency ranges. In this field, all such radiation iscolloquially referred to as THz radiation, particularly that in therange from 25 GHz to 100 THz, more particularly that in the range of 50GHz to 84 THz, especially that in the range from 100 GHz to 50 THz.

Such radiation is non-ionising and, as a result, it is particularly ofuse in medical applications. In medical imaging, the radiation isgenerally reflected from or transmitted through the patient.

Components of the sample being imaged will have a frequency dependentabsorption coefficient and refractive index, thus each component of asample subjected to radiation will leave its own characteristicfingerprint in the detected radiation. Thus, researchers have attemptedto image samples using a plurality of frequencies to create an imagefrom spectral information.

Measurements have been made using both frequency domain techniques,(where the amplitude of each frequency components is analysed), and timedomain techniques, (where the radiation is analysed as a function of thedelay time introduced by the sample into the path of the radiation).Time domain imaging is described in earlier patent GB 2 347 835.

For imaging, in the time domain, a time domain spectra is obtained foreach pixel of the sample. It is then necessary to obtain a singleparameter from each spectra to plot for each pixel. Previous attempts atproducing images from such spectra have used the amplitude of thehighest maximum or lowest minimum.

SUMMARY OF THE INVENTION

It is an object of the present invention to provide a method andapparatus which derives a parameter from time domain spectra to allow animage with enhanced contrast to be generated. The image can be a twodimensional image of an area of a sample, or it may be a profile of aline passing through the sample.

Thus, in a first aspect, the present invention provides a method ofimaging a sample, the method comprising:

(a) irradiating a sample with a pulse of electromagnetic radiation, saidpulse having a plurality of frequencies in the range from 25 GHz to 100THz;

(b) determining a first parameter at least related to the amplitude ofthe radiation, which is either reflected from and/or transmitted by thesample, in the time domain;

(c) calculating the value of the first parameter at a first time valuerelative to the value of the first parameter at a second time value,wherein the second time value coincides with a physical feature of thefirst parameter with respect to time; and

(d) generating an image by plotting the value calculated in step (c) fordifferent points of the sample.

As previously described, the image may be just a profile through thesample or a full 2D image of an area of the sample or even a 3D image ofa volume of the sample.

The inventors have found that the method of the present invention hasproduced surprisingly good results for identifying the extent of tumoursand particularly Basal Cell Carcinoma (BCC) which is the most commonform of skin cancer in caucasians. BCCs seldom metatise, but they can belocally invasive and with 80% occurring on the head and neck, theireffects on the patient can be cosmetically traumatic.

The value of the first parameter at a first time value is calculatedrelative to the value of the first parameter of the second time value.The second time value is chosen to coincide with a physical feature ofthe dataset with respect to time such as a minima, maxima, zero-crossingpoint, point of inversion, discontinuity etc.

The second time value preferably corresponds to a predetermined maximaor minima of the first parameter with respect to time.

For example, the value of the first parameter at the first time value,i.e. E(t₁) may be divided or multiplied by the value of the firstparameter at the second time value i.e. E(t₂). Alternatively, these twovalues may be added or subtracted from one another. Some furthernumerical relationship may also be implemented.

It has been found for imaging tumours and especially BCCs that goodcontrast can be obtained if Z is plotted where Z:

$Z = \frac{E\left( t_{1} \right)}{E\left( t_{2} \right)}$

Preferably, t₁ is chosen as one of the time values ascending ordescending to or from a local maxima or minima where t₂ is chosen tocoincide with this local maxima or minima.

In a particularly preferred configuration, the first time value ischosen from the subsequent time values ascending from a local minima,where t₂ is chosen to coincide with that minima.

When constructing the image, it is necessary to derive Z for each pixelof the sample of interest. t₂ is tied to a physical feature of the timedomain plot of the first parameter which may move in time from point topoint. For example, t₂ may coincide with a minima which shifts slightlyfrom point to point. The first time value t₁ may be fixed in time or maybe located at a fixed time interval δt from t₂ such that t₁ moves witht₂ from point to point of the sample.

The first time value may be chosen to coincide with the signal receivedfrom a point of interest within the sample. For example, if the sampleis measured using reflection geometry, the first value may coincide withthe time in which the reflected pulse takes the reach the part ofinterest in the sample, e.g. if the sample is a BCC sample, the timewhich it takes to be reflected from the tumour. If t₂ is chosen tocoincide with a main minima obtained using reflection geometry, δt maybe chosen as approximately the time which the radiation takes to travelfrom the sample surface to the tumour and back to the surface. Inpractice, due to variations in the composition of the sample, t₁ will bedetermined experimentally.

As the contrast can be optimised by careful selection of the first timevalue, preferably, the method further comprises the step of determiningthe first parameter for a plurality of first time values and calculatingthese values of the first parameter relative to the value of the firstparameter at the second time value. The method then further comprisingthe step of generating an image of the sample for each of the said firsttime values.

Thus, by viewing each of the images, an optimum first time value can bedetermined. Alternatively, each image may be analysed by means of acomputer in order to derive the best value of t₁ using computer means.

If t₁ is fixed for each point in a single image, then the image may beoptimised by finding the optimum value of the fixed t₁. If t₁ is alocated fixed time interval δt from t₂ for each point in the image, thenthe image may optimised by finding the optimum value of the timeinterval, δt.

Good contrast may also be obtained by careful choice of the fixed timeinterval between the first and second time values regardless of whetheror not the second time value corresponds to a physical feature of thetime domain data set. Thus, in a second aspect, the present inventionprovides a method of imaging a sample, the method comprising:

(a) irradiating a sample with a pulse of electromagnetic radiation, saidpulse having a plurality of frequencies in the range from 25 GHz to 100THz;

(b) determining a first parameter related to the amplitude of theradiation, which is either reflected from and/or transmitted by thesample, in the time domain;

(c) calculating the value of the first parameter at a first time valuerelative to the value of the first parameter at a second time value; and

(d) generating an image by plotting the value calculated in step (c) fordifferent points of the sample, wherein the time interval between thefirst and second time values is constant for all points used to generatethe image.

The value of the first parameter for the first time value with respectto the second time value may be calculated using any of the methodssuggested in relation to the first aspect.

As enhanced contrast may be obtained by measuring the time domain signalat just two points, the image acquisition time may be substantiallyreduced. However, in some cases, it may be necessary to still acquirethe data for more than two time values in order to derive the firstparameter.

The first parameter may be the amplitude of the radiation itself.However, preferably, the detected amplitude is processed in order toobtain a first parameter.

When imaging non-rigid samples, it is preferable if a flat surface isprovided for receiving the radiation. This is typically provided by amember such as a window which is transparent to the irradiatingradiation. The sample is then pressed against this window to provide aflat analysis surface. However, the window itself can cause its ownproblems. This is because the window, has two interfaces, a frontinterface upon which the irradiating radiation is directly incident anda back surface which abuts the sample. Both of these surfaces willreflect radiation given rise to at least two unwanted signals in thedetected reflected radiation. We refer to at least two reflectionsbecause multiple internal reflections may occur within the window givingeven further unwanted effects.

Preferably, if the sample is irradiated through such a member, steps aretaken in order to remove this so-called “baseline” signal. Thus, if step(a) comprises a step of irradiating the sample through a member which istransparent to the irradiating radiation and which abuts the sample, themethod preferably further comprises the step of obtaining a baselinesignal by irradiating a member in the absence of the sample anddetecting the amplitude of the reflected radiation. The first parameteris then determined by subtracting the baseline signal from the detectedamplitude of the reflected radiation from both the member and thesample.

The baseline signal can be subtracted from the amplitude of the detectedreflected radiation in the time domain or the frequency domain. Thefrequency domain signal is created by Fourier transforming the timedomain signal.

In addition to or as an alternative, it may also be desirable to obtaina reference signal which is obtained by replacing the sample with thereference object of known reflectance or transmittance and measuring theamplitude of radiation reflected from the reference object (if thesample is to be measured in reflectance mode) or the transmittance ofthe object (if the sample is to be measured in transmission mode). Thefirst parameter is then derived by dividing the detected amplitude ofthe radiation from the sample with the reference signal in the frequencydomain.

Where both a reference signal and a baseline signal are obtained, thebaseline signal is subtracted from both the sample signal and thereference signal and the baseline subtracted sample signal is thendivided by the baseline subtracted reference signal. This division isperformed in the frequency domain.

It is possible to interchange the reference and baseline signals. Forexample, the baseline signal may be obtained by replacing the samplewith an object of known reflectance or transmittance, for example, aQuartz window and the reference signal may be obtained in the absence ofa sample or a reference object.

Preferably, the amplitude of the detected radiation or the amplitude ofthe detected radiation which has been baseline subtracted and/or dividedby a reference signal is then filtered. More preferably, this isachieved by multiplying the amplitude which may have been baselinesubtracted and/or divided by a reference signal in the by the complexFourier transform of a filter function F(t), where F(t) is a non-zerofunction whose integral between time limits t_(a) and t_(b) is zero, andwhere t_(a) and t_(b) are chosen to encompass the part of interest ofthe reflected or transmitted pulse of radiation.

To clarify, radiation is supplied to the sample in the form of a pulse.The pulse will have a finite length. As the pulse passes through freespace and through the sample, the length of the pulse may change. It isimportant to set t_(a) and t_(b) so that integration is performed overthe part of the pulse of interest.

Considering the case of transmission, if the sample has the sameproperties as free space, then the radiation will pass through thesample without any delay. However, if the sample introduces a time delayinto the beam, then t_(a) will be set to a negative value which has amagnitude is equal to or larger than the expected value of this timedelay. T_(b) will generally, be set to the positive value of t_(a).

Similarly, during reflection, t_(a) will be set to a negative valuewhich has a magnitude is equal to or larger than the expected value ofthis time delay introduced by the pulse being reflected from the deepestpoint of interest in the sample. t_(b) will generally, be set to thepositive value of t_(a).

Generally, the method of the above aspect of the present invention willbe performed using the above described apparatus where the radiation ismeasured using a reference beam and wherein a scanning delay line isintroduced into the path of the reference beam or irradiating beam inorder to measure the phase change introduced by the sample. In thissituation, t_(a) and t_(b) will be the negative and positive limits ofthe scanning delay line. These may be set to the duration of the pulseor possibly a shorter time range.

Preferably function F(t) comprises a Gaussian component. Generally,t_(a) and t_(b) are symmetric about 0, preferably F(t) is also symmetricabout zero.

As the signal is usually digitally sampled, strictly a summation isperformed as opposed to an integral.

A particularly preferable form of F(t) is provided by:

${F(t)} = {\frac{2}{\pi}\left\{ {\frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\alpha})}^{2}}}{\alpha} - \frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\beta})}^{2}}}{\beta}} \right\}}$where α and β are constants.

Preferable α is substantially equal to the shortest pulse length of thebeam of pulsed radiation and β is set to be much longer than the pulselength, typically 5 to 100 times the pulse length. However, both ofthese values will generally be optimised by the operator.

If β is greater than or comparable to the time which it takes theradiation to penetrate the sample to the point of interest then, F(t)can take the simplified form:

${F(t)} = {\frac{2}{\pi}\left\{ {\frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\alpha})}^{2}}}{\alpha} - \frac{1}{T}} \right\}}$where α is a constant which is substantially equal to the shortest pulselength of the beam and T is substantially equal to the time which ittakes the beam of radiation to penetrate to the deepest point ofinterest in the sample.

A sample signal which has been baseline subtracted, divided by areference signal and filtered is referred to as the impulse function.Preferably, the first parameter is the impulse function.

As the present invention is concerned with producing an image, theactual value of the impulse function at a pixel does not have to have astrong numerical correlation with a physical parameter of the sample.However, it is important that the impulse function is capable of showingvariations in the composition of the sample.

The first time value may be determined by:

(i) selecting n regions of said sample, where n is an integer of atleast 2;

(ii) producing n time domain spectra of a second parameter which is atleast related to the amplitude of the radiation reflected from and/ortransmitted by the sample, where the nth spectra is obtained byaveraging the second parameter across the nth area and plotting thisaveraged value, as it varies with delay time;

(iii) determining the first time value by comparing at least two of then spectra derived in step (ii).

Preferably, each spectra is normalised to the value of the secondparameter at a time value which coincides with a physical feature of thesecond parameter with respect to time, prior to comparing the spectra. Aphysical feature of the spectra of the second parameter may be aminimum, maximum, point of inflection etc. Preferably, the same physicalfeature is chosen which identifies the second time value.

The spectra may be compared by performing a mathematical operation on atleast two of the spectra. For example, the spectra may be divided by orsubtracted from one another.

The selected regions may be lines or areas of the sample. They may becontinuous lines or areas or discontinuous lines or areas. The regionsmay be of the same size of different sizes and may comprise any numberof pixels.

Preferably, the regions selected by first generating a preliminary imageof an area of the sample by plotting a preliminary imaging parameterwhich is at least related to the amplitude of the reflection reflectedfrom and/or transmitted by the sample.

The preliminary imaging parameter may be any parameter which can be usedto generate an image of the sample, for example, it may be the value ofthe first parameter and/or a second parameter, at a minima, maxima orpoint of inflection, etc. Alternatively, it may be the value of thefirst parameter or second parameter at an arbitrary time value or a timevalue which has previously been used to show contrast in the sample orsimilar samples.

In the simplest case, just two areas are selected. However, the methodmay also be used to optimise contrast between 3 or more differentregions. To achieve this, three regions are chosen, the value of thepreliminary imaging parameter differing between said three selectedregions, the first time value being chosen to coincide with a regionwhere there is a difference between the three spectra.

The first parameter and second parameter may be the same or different.The second parameter may be preferably chosen from any of the parameterspreviously described in relation to the first parameter.

Regions of the sample may be chosen at random or, alternatively, regionsof the sample may be chosen where some contrast has already beenobserved. In this example, the n areas are chosen such that the level ofthe preliminary imaging parameter within each region is substantiallyconstant and the level of the preliminary imaging parameter is differentbetween at least two of the regions.

Preferably, a standard deviation of about 30% in the primary imagingparameter is allowed within a selected region, more preferably 20%, evenmore preferably 10%.

The areas may be selected randomly by a person. Alternatively, the areaselection may be automated. For example, a processor may be used toselect n regions at random or sequentially and then compare the spectrafrom these regions in order to optimise the time value. The processor ispreferably configured to select each area such that the preliminaryimaging parameter is substantially constant within each area aspreviously described. The processor may further be configured to varythe size of the selected region in order to select regions which fallwithin the above criteria.

The above method of deriving a first time value may also beadvantageously used with the method described in GB 2 360 186. Thus, ina third aspect, the present invention provides a method of imaging asample, the method comprising:

(a) irradiating a sample with a pulse of electromagnetic radiation, saidpulse having a plurality of frequencies in the range from 25 GHz to 100THz;

(b) determining a first parameter at least related to the amplitude ofthe radiation, which is either reflected from and/or transmitted by thesample, in the time domain; and

(c) generating an image of the sample by plotting the value of the firstparameter at a first time value, said first time value being determinedby:

-   -   (i) selecting n regions of said sample, where n is an integer of        at least 2,    -   (ii) producing n time domain spectra of a second parameter which        is at least related to the amplitude of the radiation reflected        from and/or transmitted by the sample, where the nth spectra is        obtained by averaging the second parameter across the nth area        and plotting this averaged value, as it varies with delay time;    -   (iii) determining the first time value by comparing at least two        of the n spectra derived in step (ii).

In a fourth aspect, the present invention provides imaging apparatuscomprising:

a source for irradiating a sample with a pulse of electromagneticradiation, said pulse having a plurality of frequencies in the rangefrom 25 GHz to 100 THz;

a detector for detecting the amplitude of radiation which is eitherreflected from or transmitted by the sample;

means for determining a first parameter which is related to theamplitude of the reflected and/or transmitted radiation in the timedomain;

calculating means for calculating the value of the first parameter at afirst time value relative to the value of the first parameter at asecond time value, wherein the second time value coincides with aphysical feature of the dataset of the first parameter with respect totime; and

means for generating an image by plotting the value calculated by thecalculating means for different points of the sample.

In a fifth aspect, the present invention provides an apparatus forimaging a sample, the apparatus comprising:

a source for irradiating a sample with a pulse of electromagneticradiation, said pulse having a plurality of frequencies in the rangefrom 25 GHz to 100 THz;

detector means for detecting the amplitude of radiation reflected fromand/or transmitted by the sample;

means for determining a first parameter related to the amplitude of theradiation in the time domain;

calculating means for calculating the value of the first parameter at afirst time value with respect to the value of the first parameter at asecond time value, and

imaging means for generating an image by plotting the value calculatedby the calculating means for different parts of the sample, wherein thetime interval between the first and second time intervals is constantfor all points used to generate the image.

The detector may be a direct detector of THz radiation or an apparatusaccording to any of the above three aspects of the invention, thedetector may be a direct detector of THz radiation or it may be of thetype which converts THz radiation into an easily readable signal.

For example, the detector may comprise a non-linear crystal which isconfigured such that upon irradiation of a probe beam and a THz beam,the polarisation of the probe beam is rotated. The probe beam can be ofa frequency which can be easily measured (for example near infra-red).Typical crystals which exhibit this effect, the so-called “AC Pockels”effect are GaAs, GaSe, NH₄H₂PO₄, ADP, KH₂PO₄, KH₂ASO₄, Quartz, AlPO₄,ZnO, CdS, GaP, BaTiO₃, LiTaO₃, LiNbO₃, Te, Se, ZnTe, ZnSe, Ba₂NaNb₅O₁₅,AgAsS₃, proustite, CdSe, CdGeAs₂, AgGaSe₂, AgSbS₃, ZnS, organic crystalssuch as DAST (4-N-methylstilbazolium. This type of detection mechanismis generally referred to as ‘Electro-optic sampling’ or EOS.

Alternatively, the detector could be a so-called photoconductingdetector. Here, the detector comprises a photoconductive material suchas low temperature grown GaAs, Arsenic implanted GaAs or radiationdamaged Si on Sapphire. A pair of electrodes, for example in a bow-tieconfiguration or in a transmission line configuration are provided on asurface of the photoconductive material. When the photoconductivematerial is irradiated by the reflected radiation and also, the probebeam, a current is generated between the two electrodes. The magnitudeof this photovoltage current is an indication of the magnitude of theTHz signal.

Although it is possible to generate THz radiation directly, the mosteffective THz generation can be achieved by converting a pump beam intoa THz beam. To do this, the source comprises a frequency conversionmember and a source of a pump beam.

There are many possible options for the frequency conversion member. Forexample, the frequency conversion member may comprise a non-linearmember, which is configured to emit a beam of emitted radiation inresponse to irradiation by a pump beam. Preferably, the pump beamcomprises at least two frequency components, (or two pump beams havingdifferent frequencies are used), the non-linear member can be configuredto emit an emitted beam having a frequency which is the difference ofthe at least two frequencies of the pump beam or beams. Typicalnon-linear members are: GaAs or Si based semiconductors. Morepreferably, a crystalline structure is used. The following are furtherexamples of possible materials:

-   NH₄H₂PO₄, ADP, KH₂PO₄, KH₂ASO₄, Quartz, AlPO₄, ZnO, CdS, GaP,    BaTiO₃, LiTaO₃, LiNbO₃, Te, Se, ZnTe, ZnSe, Ba₂NaNb₅O₁₅, AgAsS₃,    proustite, CdSe, CdGeAs₂, AgGaSe₂, AgSbS₃, ZnS, GaSe or organic    crystals such as DAST (4-N-methylstilbazolium).

In order to produce an emitted beam having a frequency in the THzregime, preferably the at least two frequencies of the pump beam orbeams are in the near infra-red regime. Typically, frequencies between0.1×10¹² Hz and 5×10¹⁴ Hz are used.

Alternatively the frequency conversion member is a photoconductingemitter, such an emitter comprises a photoconductive material such aslow temperature grown or arsenic implanted GaAs or radiation damaged Sior Sapphire.

Electrodes which may be of any shape such as a dipole arrangement, adouble dipole arrangement, a bow-tie arrangement or transmission linearrangement are provided on the surface of the photoconductive material.At least two electrodes are provided. Upon application of a bias betweenthe electrodes and irradiation of a pump beam(s) having at least twodifferent frequency components, a beam of radiation is emitted having afrequency different to that of the at least two frequency components ofthe pump beam or beams.

When a pulse having a plurality of frequencies passes via a sample to adetector, the various frequencies will not arrive at the detector at thesame time due to the frequency dependent response of the sample. A timedomain signal can be established by measuring the amplitude of thedetected radiation with respect to time. In order to achieve this, it ispreferable if a scanning delay line is inserted into either the path ofthe probe or pump beam. The delay line can be configured to scan overthe whole length of the pulse.

As mentioned with respect to the method of the present invention, thefirst time value can be selected in order to optimise contrast in theimage. The first time value may be fixed with respect to zero across animage, or it may be fixed with respect to the second time value.

Therefore, preferably, the apparatus further comprises means forgenerating a plurality of images corresponding to different first timevalues, the apparatus also comprising means for scanning through saidimages by varying said first time value in order to determine a firsttime value which allows the best contrast.

The apparatus preferably further comprises:

(i) means for selecting n regions of said sample, where n is an integerof at least 2;

(ii) means for producing n time domain spectra of a second parameterwhich is at least related to the amplitude of the radiation reflectedfrom and/or transmitted by the sample, where the nth spectra is obtainedby averaging the second parameter across the nth area and plotting thisaveraged value as it varies with delay time; and

(iii) means determining the first time value by comparing at least twoof the n spectra.

In a sixth aspect, the present invention provides:

(a) means for irradiating a sample with a pulse of electromagneticradiation, said pulse having a plurality of frequencies in the rangefrom 25 GHz to 100 THz;

(b) means for determining a first parameter, at least related to theamplitude of the radiation, which is either reflected from and/ortransmitted by the sample, in the time domain; and

(c) means for generating an image of the sample by using the value ofthe first parameter at a first time value,

said apparatus being configured to determine said first time value by:

-   -   (i) selecting n regions of said sample, where n is an integer of        at least 2;    -   (ii) producing n spectra in the time domain, where the nth        spectra is obtained by plotting the second parameter, averaged        across the nth area, against delay time;    -   (iii) determining the first time value by comparing at least two        of the n spectra.

The means for selecting regions may be automatic means. For example, aprogram may be executed to select two or more regions of the sample atrandom or two or more regions of the sample which have the largestdifference between their mean colour levels, the standard deviation ofthese areas may be calculated to ensure that the selected area does notcover an area of high contrast.

Alternatively, the apparatus may be provided with means which allow anoperator to select two or more areas. For example, the apparatus may beconfigured to support a cursor which the operator may use to selectparticular areas of the image.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention will now be described with reference to thepreferred non-limiting embodiments in which:

FIG. 1 is a schematic of a reflection imaging system which may be usedin accordance with an embodiment of the present invention;

FIG. 2 is a is a schematic of a transmission imaging system which may beused in accordance with an embodiment of the present invention;

FIG. 3 a is a perspective view of a sample holder and FIG. 3 b is a planview of a sample holder which may be used in accordance with anembodiment of the present invention;

FIG. 4 illustrates a plurality of deconvoled impulse functions for asample of a skin tumour;

FIG. 5 is a schematic illustrating the analysis of an impulse functionin accordance with an embodiment of the present invention;

FIG. 6 a illustrates a visible image of a skin tumour, FIG. 6 billustrates a corresponding Terahertz image of the same sample producedin accordance with an embodiment of the present invention and FIG. 6 cis a plot of the mean value of the THz signal of the image of FIG. 6 bplotted for four pixels of the image;

FIG. 7 a is a vertical histology section of a skin tumour and FIG. 7 bis a corresponding section through a Terahertz image obtained inaccordance with an embodiment of the present invention;

FIG. 8 is a schematic illustrating the analysis of an impulse functionin accordance with a further embodiment of the present invention;

FIG. 9 is a schematic illustrating the analysis of an impulse functionin accordance with a yet further embodiment of the present invention;

FIG. 10 is a schematic illustrating the analysis of an impulse functionin accordance with an even further embodiment of the present invention;

FIG. 11 is an apparatus which can optimise the contrast of a THz imagein accordance with an embodiment of the present invention;

FIG. 12 is an image of a tumour generated using E(t)/E(min) as explainedwith reference to FIG. 5;

FIG. 13 illustrates two time domain spectra averaged over areas “a” and“b” of FIG. 12;

FIG. 14 illustrates a spectra derived from the difference between thetwo spectra of FIG. 13;

FIG. 15 illustrates a spectra derived by dividing one spectra of FIG. 13by the other spectra of FIG. 13;

FIG. 16 illustrates a plot of two time domain spectra derived byplotting E(t)/E(min), the spectra are similar to those of FIG. 13 buthave been normalised to the value E(min);

FIG. 17 illustrates a time domain spectra derived from the difference ofthe two spectra of FIG. 16;

FIG. 18 illustrates a time domain spectra derived from dividing onespectra of FIG. 16 by the other spectra of FIG. 16;

FIG. 19 illustrates an image of a tumour derived by plotting themagnitude of the main minima of the impulse function for each pixel inthe x and y direction;

FIG. 20 illustrates the mean value of the minimum view averaged overeach of the areas s1, s2, d1, d2, h1 and h2 indicated in FIG. 19;

FIG. 21 illustrates a plot of the mean impulse function in the timedomain for areas s1, s2, d1, d2, h1 and h2 of FIG. 19;

FIG. 22 illustrates a plot similar to that of FIG. 21 where the spectrahave been normalised to the value E(min);

FIG. 23 illustrates a detail of the plot of FIG. 21;

FIG. 24 illustrates a detail of the plot of FIG. 22;

FIG. 25 is a plot identical to that of FIG. 21 but repeated to alloweasy comparison;

FIG. 26 is a plot identical to that of FIG. 22 but repeated to alloweasy comparison;

FIG. 27 illustrates a detail of the plot of FIG. 25;

FIG. 28 illustrates a detail of the plot of FIG. 26;

FIG. 29 illustrates a spectra of two time domain spectra, the upperspectra being produced by subtracting the spectra obtained from area h1in FIG. 19 from the spectra obtained from the diseased tissue areas d1and d2; the lower plot is a spectra obtained by subtracting a meanspectra obtained from healthy tissue area h1 of FIG. 19 from inflamedtissue areas s1 and s2;

FIG. 30 illustrates an enhanced image obtained by plotting the parameterE(t)/E(min) for t=1.2 picoseconds; and

FIG. 31 illustrates an image obtained by plotting E(t)/E(min) for t=3.7picoseconds.

FIG. 1 shows a basic THz reflection imaging system. The system can besplit into three main sections: a generation section 1, an imagingsection 3 and a detection section 5.

A THz beam of radiation is produced in generating section 1. Generationsection 1 comprises an optically amplified Ti:Sapphire laser 7 toproduce a beam 9 of 250 fs pulses at 800 nm.

The optically amplified laser 7 comprises a solid state pump laser 11and an amplified Ti:Sapphire section 13. Beam 9 is divided by 50:50 beamsplitter 15. Beam splitter 15 divides the beam into a pump pulse 17which is used to generate THz radiation and a probe pulse 19 which isused in the detection of the radiation reflected from the sample.

THz pulses are generated by optical excitation and charge accelerationof biased GaAs stripline antenna 21 which outputs THz pulses 23 whichhave a bandwidth of 0.1 THz to 3 THz. The emitted beam 23 is then passedinto imaging section 3.

The irradiating THz beam 23 is then directed via first off-axisparabolic mirror 25 onto second off-axis parabolic mirror 27 and thenfocussed to a 400 μm spot on sample 29. Sample 29 is mounted onmotorised stage 31. The sample mount 31 is computer controlled and canbe configured to move the sample such that each section (pixel) of thesample in question can be imaged.

The reflected radiation is then collected via third off-axis parabolicmirror 33 and reflected via planar mirror 35 into the detection system5. In the detection system, the reflected radiation 37 is reflectedoff-axis parabolic mirror 39. Off-axis parabolic mirror 39 has a holethrough it which allows the reflected beam 37 to be combined with theprobe beam 19 obtained by beam splitter 15.

Prior to combining the probe beam 19 with the reflected radiation 37 atoff-axis parabolic mirror 39, the probe beam 19 is reflected off planarmirror 41 into scanning delay line 43. Scanning delay line comprisescuboid mirror 45 which moves back and forth under computer control 47 inorder to vary the pulse length of the probe beam 19 with respect to thepump beam 17 and hence the irradiating and reflected radiation. Thus,the relative phase of the probe beam 19 can be varied with respect tothe pump pulse 19.

The combined probe beam and reflected radiation 37 is then collimatedand focused onto ZnTe detection crystal 49. The THz radiation isdetected using the linear electro-optic Pockel's effect.

The linearly polarised probe beam 19 has its polarisation oriented suchthat it has components along both the ordinary and extra-ordinary axisof detection crystal 49. Each of the axis has distinct refractiveindices n₀ and n_(t) along the ordinary and extra-ordinary axis of thecrystal 49 respectively. In the absence of reflected radiation 37, thelinearly polarised probe beam 19 passes through the detection crystalwith a negligible change in its polarisation.

The applicant wishes to clarify that the angle through which thepolarisation is rotated is negligible. However, the linearly polarisedbeam can become slightly elliptical. This effect is compensated for by avariable retardation waveplate, for example, quarter waveplate 51.

The beam emitted by detection crystal 49 is converted into circularlypolarised beam by quarter waveplate 51. The beam is then split into twolinearly polarised beams by Wollaston prism 53 (or an equivalent devicefor separating orthogonal polarisation components) which directs the twoorthogonal components of the polarised beam onto a balance photo diodeassembly 55. The balance photo diode signal is adjusted using waveplate51 so that the difference in outputs between the two diodes is zero whenonly the probe beam 19 is incident on the detection crystal 49.

However, if the detector also detects reflected radiation 37, the anglethrough which the polarisation is rotated by is not negligible. This isbecause the THz electric field modifies the refractive index of theprobe beam 19 along one of the axis n_(e), n_(o). This results in thephysical field emitted by detection crystal 49 being elliptical.Therefore, the polarisation components separated by prism 53 are notequal and the difference in voltage between the output diodes of balancephoto-diode assembly 55 gives a detection voltage.

The rotation in the polarisation is obscured if the phase of thereflected radiation 37 reaching the detection crystal 49 is different tothe phase of the probe beam 19. As the sample 29 will introduce adifferent phase change for different frequency components of theirradiating radiation 23, the delay line 43 is swept under control ofcomputer 47 in order to sweep the phase of the probe beam 19 withrespect to the reflected radiation 37 and hence detect the delay time ofeach of the frequency components of reflected radiation 37. Thisproduces a so-called time-domain signal.

Computer control 47 also is used to control the movement of stage 31 inorder to ensure that a full data set is collected for each pixel ofinterest of sample 29.

The above system illustrates is configured for reflection measurements.However, it will be appreciated that those skilled in the art that thesystem could also be adapted for a transmission method. Such anapparatus is shown in FIG. 2. To avoid unnecessary repetition, likereference numerals will be used to denote like features. In imagingsection 3, first 25 and second 27 parabolic mirrors are configured tofocus radiation onto sample 29. The third parabolic mirror 23 thencollects the transmitted radiation and directs it into detection section5. The generation 1 and detection 5 sections are the same as those forthe reflection apparatus of FIG. 1.

The method of the present invention has been found to be of particularuse in imaging Basal Cell Carcinoma. FIG. 3 illustrates a preferredmounting method for a sample of this type. FIG. 3 a shows a perspectiveview of the sample holder whereas FIG. 3 b shows a cross section of thesample holder.

FIGS. 3 a and 3 b show two samples 61 and 63 located in a hollow,cylindrical sample holder 65. The sample holder 65 is based on amodified liquid cell of the type used for Fourier transform infra-redspectrometry. The liquid cell is cylindrical in shape and comprises thehollow cylindrical container 65, a 2 mm thick Quartz plate 67 with adiameter of 25 mm, is placed at the base of container 65. The THzradiation irradiates the sample through Quartz window 67.

As water is a strong absorber of radiation in the THz regime, excessfluid is removed from the sample using lint free paper. The samples weremounted on the Quartz window 67 using tweezers and the top surface ofthe skin was placed on Quartz window 67. Sponge 69 is then placed on topof the sample so that the samples are sandwiched between the sponge 69and Quartz window 67. The retaining ring 71 is placed on top of Quartzwindow 67 in order to keep the window level and also to maintain theposition of the samples 61 and 63.

A polythene window 73 is provided on the opposing side of the sponge tothe Quartz window and in contact with the sponge. The polythene windowis held in place by retainer plate 75 which is screwed into the sampleholder 65 by screw 77. This construction provides an airtight holder forthe sample which prevents the sample from drying out during imaging.

FIG. 4 illustrates a time domain impulse function for diseased andhealthy tissue. FIG. 4 shows the impulsed function along the y axisplotted against time delay (arbitrary units) along the x axis. Fivetraces are shown in FIG. 4 for five different pixels.

For each pixel, the entire THz waveform is recorded and averaged overfour readings to reduce any fluctuations arising from the laser. Thesignal to noise ratio is approximately 1,000 to 1.

The impulse function is derived from the amplitude of the reflectedradiation. In order to derive the input function for a particular pixel,two further measurements are made in addition to a measurement of thesample.

First, a baseline signal is measured. In the reflection measurement, itis important to try to reduce the effects of spurious reflections. Theapparatus of FIG. 3 uses a Quartz window 67 through which the samples 61and 63 are illuminated. Quartz window 67 has a lower interface 79 whichwill cause spurious reflections and also an upper interface 81 whichwill cause spurious reflections.

It is important to try to eliminate the effects of these two interfacesfrom the measurements. Therefore, prior to imaging the sample, a secondQuartz window (not shown) is placed in the sample holder instead ofsample 61 and 63. The amplitude of the detected THz signal is measuredfor each pixel so that a baseline signal is obtained through each pixel.Variations in Quartz window 67 and also variations in the incident angleof the irradiating radiation may cause artificial variations of themeasurements of the sample. Therefore, the baseline is obtained for eachpixel.

The baseline measurement was averaged twenty times to obtain baselinesignal B(t).

A reference spectrum is also measured. This is measured by using thesample bolder with just Quartz window 67 on its own. This referencesignal will be represented in the time domain as R(t).

First, in order to obtain an impulse function for one waveform, thebaseline signal B(t) is subtracted from the measured time domainwaveform S(t). This operation can be formed in either the time domain orfrequency domain.

This baseline subtracted waveform should provide data where the spuriousreflections due to the lower 79 and upper 81 interfaces of Quartz window67 have been eliminated.

Next, the baseline subtracted waveform is complex Fourier transformed togive:S′(v)−B′(v)

In order to use the reference signal, the baseline signal is againsubtracted from the reference waveform R(t). This subtraction operationcan be performed in either the time domain or the frequency domain. Theresult is then complex Fourier transformed to give:R′(v)−B′(v)

The baseline subtracted measured data is then divided by the baselinesubtracted reference signal in the frequency domain to give:

$\frac{{S^{\prime}(v)} - {B^{\prime}(v)}}{{R^{\prime}(v)} - {B^{\prime}(v)}}$

The data is then filtered. In general, a filter function will bedescribed by F(t) and is complex Fourier transform will be representedby a F′(v).

A filter function preferably is used because the THz pulse system cangenerate and detect pulses comprising frequencies over some finiterange, typically from less than 100 GHz to over 3 Thz. There is a highfrequency limit above which the THz signal falls below the noise levelof the detection system.

Similarly, the THz signal level falls below the noise level at lowfrequencies. Thus, there is a need to remove the high and low frequencynoise. A particularly preferable function for achieving this is:

${F(t)} = {❘{\frac{2}{\pi}\left\{ {\frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\alpha})}^{2}}}{\alpha} - \frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\beta})}^{2}}}{\beta}} \right\}}}$

The parameters α and β are selected to control the high and lowfrequency roll-off of the function. α is set approximately the shortestTHz pulse length (half cycle) obtainable within the THz system. β is setto be much longer than the THz pulse. In operation, the two parametersare optimised manually by the operator to obtain the best compromisebetween bandwidth and noise.

As the above function comprises two Gaussian functions with similarareas but opposite signs, the above function ensures that the integralof the filter function for all time is zero.

If the value of β is comparable to or greater than the total time-delayscan range, then an alternative function can be used:

${F(t)} = {❘{\frac{2}{\pi}\left\{ {\frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\alpha})}^{2}}}{\alpha} - \frac{1}{T}} \right\}}}$where T represents the total range of delay times used i.e.T=T_(max)−T_(min). This ensures that the overall integral from T_(min)to T_(max) is always zero.

The filter function is multiplied by the complex Fourier transform ofthe signal which is being corrected by the baseline and the reference.

Once the impulse function has been derived as illustrated in FIG. 4, theimage is then generated using the analysis described with reference toFIG. 5.

FIG. 5 shows a schematic impulse function for just one pixel. In orderto generate an image of the sample, it is important to derive ameaningful parameter from this trace which can be plotted for eachpixel. By a meaningful parameter, we mean a parameter which can be usedto indicate and emphasise variations in the sample.

The inventors have surprisingly found that the following method ofderiving a parameter provides an image with a remarkable sharp contrastfor illustrating a number of variations in the sample but has found tobe of particular use for indicating the presence of Basal CellCarcinoma.

The value of the impulse function at the minima is measured. Also, thevalue of the impulse function at a time t. In this example, the time tis chosen from the times which correspond to ascending values of thefirst parameter after the minima.

A parameter Z is derived by dividing the first parameter at a time t bythe value of the first parameter at the minima i.e.

$Z = \frac{E(t)}{E\left( \min \right)}$

To generate an image, each value of Z is assigned a colour and a fullcolour map is generated by plotting z on the z axis against x and ywhich define the pixel positions. As the contrast varies with the chosenvalue of t, t is chosen in order to give the best contrast. “t” is fixedfor the image.

FIG. 6 a illustrates a visual image of a non-recurrent tumour present ona nasal tip. The tumour is 7 mm in diameter. The tumour 91 is shown onthe left of the image with healthy tissue 93 shown on the front of theimage. A histology cut is indicated on both diagrams 6 a and 6 b with adotted line. The importance of this will be described with reference toFIG. 7.

The false colour THz image 6 b shows a strong contrast between thediseased and healthy tissue. The healthy tissue 93 shows a fairlyuniform greyscale contrast throughout. The strong contrast between thediseased and healthy tissue which is illustrated in FIG. 6 b is notpresent in the visible image where it is difficult to differentiatebetween the two.

Four areas, d1, d2, h1 and h2 corresponding to two areas of diseasedtissue and two areas of healthy tissue are indicated on FIG. 6 b. FIG. 6c illustrates a plot of the main value of the impulse function for eachof these areas. It can be seen that the mean value for the two diseasedareas is much higher than that of the mean value for the two healthyareas. The standard derivation is used to illustrate an error bar foreach of the regions. It can be seen that even in the worse casescenario, it is possible to clearly distinguish the diseased areas d1and d2 from healthy areas h1 and h2.

A histological analysis of the sample was then performed. A verticalhistology section was taken through the sample illustrated in FIGS. 6 aand 6 b along the dotted line. The results of which are indicated inFIGS. 7 a and 7 b.

In the histology section shown in FIG. 7 a, the epidermis is shown as athick line 101 of histology section 7 a. It is noted that the epidermisis missing in the centre of the section. This is because the tissue ispreferentially scraped prior to excision. As a result of increasedvascularisation around the tumour site and the changing composition ofthe skin structure within this region, the epidermis is more delicateand susceptible to preferential bleeding in the surrounding healthytissue. This effect is utilised by the surgeons in assisting identifyingthe location of the tumour. The tumour had not broken through to thesurface before preferential scrapping.

Tumour shrinkage is due to desification, a result of the histologicalpreparation process. The tumour is highly cellular and has a fibrouscasing. The epidermis is pushed to one side by the tumour and as aresult, thickened.

Turning to the THz image of FIG. 7 b, the upper line indicated on theplot is the mean value for the diseased tissue. The lower line is themean value for the healthy tissue. The THz image is taken in the samedirection as that of the histology section. In other words, it is just aslice through the greyscale image of FIG. 6 b taken along the dottedline. Towards the right hand side of the trace, diseased tissue is seenwhich corresponds to the histology section. Towards the left hand sideof the trace, the healthy tissue is seen as shown in the histologysection.

In the analysis technique shown in FIG. 8, the impulse function isderived in the same manner as that described with reference to FIG. 5.The parameter Z is then derived by subtracting the value of the firstparameter at the minima from the value of the first parameter at anarbitrary time t which occurs at a time δt after the minima i.e.Z=E(min+δt)−E(min).

The minima may move in time from point to point across the sample. t isdefined relative to E(min), thus t moves with the minima. To optimisethe image, the parameter δt is optimised as δt is fixed for each image.

FIG. 9 shows a variation on the analysis technique described withreference to FIG. 8. Again, the impulse function is derived in the samemanner as described with reference to FIG. 5.

The parameter Z is calculated by subtracting the value of the firstparameter at a time t₂ from the value of the first parameter at time t₁where t₁=t₂+δt. t₂ does not have to coincide with the minima. As shownin FIG. 9, t₁ and t₂ may lie on opposing sides of the minima or they maylie on the same side.

Again, as described with reference to FIG. 8, since t₁ is definedrelative to t₂, δt is optimised.

FIG. 10 shows a yet further technique for deriving Z. Again, asdescribed with reference to FIG. 5, the impulse function is derived as afunction of time. Z is then calculated by dividing the value of thefirst parameter at a time t which equals t_(min)+δt divided by the valueof the first parameter at the minima. In other words:

$Z = \frac{E\left( {\min + {\delta\; t}} \right)}{E\left( \min \right)}$

This technique is similar to that described in relation to FIG. 5.However, in the analysis technique of FIG. 5, t₁ is fixed across asingle image. If the minima moves (and hence t₂ moves from pixel topixel) t₁ still remains fixed. In the analysis technique described withreference to FIG. 10, t₁ is defined in relation to t₂, thus t₁ moveswith t₂ and δt is fixed for a particular image.

As mentioned with reference to FIGS. 5, 8, 9 and 10, careful selectionof t enhance the contrast. Generally, the operator of such equipmentwill be able to scan through a number of images at different t in orderto obtain the best contrast. Such an apparatus for performing this isillustrated in FIG. 11.

FIG. 11 shows a schematic view of such a system. The system comprises acontroller 111, a THz imaging system 113. The THz imaging system is ofthe type described with reference to either of FIGS. 1 and 2. To avoidunnecessary repetition, details of it will not be repeated here and thewhole system will be represented by box 113. Controller 111 has a screen115 which the operator can use to view a false colour THz image of thesample under test as explained with reference to FIGS. 4 and 5 and asillustrated in FIG. 6 b. The controller 111 has two modes of operation.A set-up mode and a fixed imaging mode. The controller switches betweenthese two modes using switch 117.

In the set-up mode, the full THz spectrum and hence impulse function isderived for each pixel of the sample. The image is then derived bycalculating Z for a given t. t may be a fixed value for all pixels ormay be fixed from a second time value which may correspond to a physicalfeature of the data set. The user can then turn dialler 119 and watchthe contrast of the image displayed on screen 115 change as t or δt isvaried. The operator can thus choose an optimum value of t or δt.

Alternatively, the controller itself may be provided with its own meansfor determining the optimum contrast as a function of t. For example,the controller can calculate the mean value of Z for each pixel andchoose a value of t where the variation in Z is at its largest.Alternatively, more complex statistical techniques may be employed whereonly representative values of the impulse function are considered indetermining the best contrast.

Once the operator or controller itself 111 has determined the bestcontrast, the operator can then switch 117 to switch the controller 111onto fixed mode. In this fixed mode, the controller will instructimaging apparatus 113 to only take data at the minima and at thepreferred value of t. This technique substantially reduces theacquisition time for the image. Once t has been determined for a sample,it can be used for that whole sample. It has also been found that thecell value of t can be used as samples of the same type. For example,for samples which are Basal Cell Carcinoma, the operator would expect touse the same value of t.

Controller 111 can be pre-programmed with t values which are known toprovide optimum contrast for certain samples. Alternatively, controller111 may be provided with a look-up table such that the operator canenter the type of sample which is to imaged and the controller canautomatically select the correct t or δt value.

FIG. 12 illustrates an image of a skin tumour obtained by dividing themagnitude of the impulse function at a time t=1098 (arbitrary units)with the impulse function at the minima of the electric field. This isexplained in detail with reference to FIG. 5. Specifically, the value ofE(t)/E(min) for t=1098 is plotted for each pixel to produce the image.

This figure will now be used to explain how the image may be optimised.Two areas “a” and “b” are selected from the image. Area “a” is believedto represent diseased tissue whereas area “b” is believed to representhealthy tissue.

The mean impulse function for both areas “a” and “b” is plotted for eachtime point as shown in FIG. 13. Formally, for area “a”, the function:Ē_(a)(t)is plotted where

${\overset{\_}{E_{a}}(t)} = \frac{\sum\limits_{i}^{n}{E_{i}(t)}}{n}$

E being the impulse function and n being the number of pixels in area“a”. The mean impulse function for area “b” is calculated in the samemanner and plotted.

The width of the curves represents the standard deviation since thewaveform is calculated over an area.

These two waveforms are then compared in order to determine the timevalues where a plot of E(t) will produce an image with the best contrastbetween areas “a” and “b”.

FIG. 14 illustrates the case where the spectra for the diseased tissueis subtracted from the waveform for the healthy tissue. It can be seenthat the largest difference between these two waveforms occurs close tot=1100. Thus, indicating that an image produced using this value of “t”as the first time value will produce good contrast. For example, theimage may be produced by plotting E(t), E(t)/E(min), E(t)−E(min) etc.

FIG. 15 illustrates the case where the two waveforms are compared bydividing the waveform for healthy tissue with that from diseased tissue.This comparison method also allows a value of “t” to be determined whichcan be used as the first time value.

In FIG. 13 the mean impulse function was derived for two areas and usedto determine an optimum value of “t”. In FIG. 16, the mean impulsefunctions of FIG. 13 are normalised to the value E_(min), formally:

$\frac{\overset{\_}{E_{a}}(t)}{E_{\min}}$is plotted, where:

${{\overset{\_}{E}}_{a}(t)} = \frac{\sum\limits_{i}^{n}{E_{i}(t)}}{n}$

E being the impulse function and n being the number of pixels in area a.The normalised mean impulse function for area 'b is calculated in thesame manner and plotted.

FIG. 17 illustrates the spectra of FIG. 16 subtracted from one anotherwhereas FIG. 18 illustrates the spectra from FIG. 16 divided by oneanother.

In the same manner as described with reference to FIGS. 14 and 15,values of “t” may be derived from the spectra of FIGS. 17 and 18 whichproduce an image with good contrast.

The present invention may also be used to compare three different typesof tissue within the sample and different t values may be selected toidentify different types of tissue.

In FIG. 19, six areas are selected, h1 and h2 which represent healthytissue, s1 and s2 which represent inflamed tissue and d1 and d2 whichrepresent diseased tissue.

The image of FIG. 19 is generated by plotting the minima of the impulsefunction for each pixel in the x and y directions i.e. Emin is plotted.The minimum mean impulse function for each of these areas, averagedacross each area, is then plotted and the graph shown in FIG. 20. Thevertical error bars represent the standard deviation.

FIG. 21 illustrates the mean impulse function E(t) for healthy tissue(calculated from areas h1, h2), inflamed tissue (calculated from areass1 and s2) and diseased tissue (calculated from areas d1 and d2) derivedfrom the sample illustrated in FIG. 19. FIG. 22 illustrates the spectraof FIG. 21 which has been normalised to the value of E_(min) i.eE(t)/E(min). The width of the spectra represents the standard deviation.

It can be seen by viewing FIGS. 20 and 21 side by side that the shape ofthe waveforms change between a plot of mean E(t) and mean E(t)/E_(min).The changes between these two waveforms may be used to distinguishbetween different tissue types.

FIGS. 23 and 24, illustrate in detail the region between t=0.8 ps andt=1.6 ps for FIGS. 21 and 22 respectively. In FIG. 23, it can be seenthat there is little difference between the curves for the diseased andinflamed tissue which just start to deviate from one another at aboutt=1.1 ps. However, the spectra which represents healthy tissue isclearly separated from the spectra from both diseased and inflamedtissue over the whole of the selected time range.

In FIG. 24, the spectra from the healthy and inflamed tissue is seen tofollow a similar trajectory whereas the spectra from the diseased tissuefollows a lower path which is well separated from the other two spectra.

The time value of 1.2 ps is clearly marked on both FIGS. 23 and 24.Generating an image by plotting E(t=1.2 ps) can thus be used to clearlydistinguish between healthy tissue on the one hand and diseased andinflamed tissue on the other. Generating an image by plotting E(t=1.2ps)/E_(min) can be used to clearly distinguish between diseased tissueon the one hand and health and inflamed tissue on the other.

By choosing a different time value, it is possible to distinguishbetween diseased and inflamed tissue. FIGS. 25 and 26 are identical toFIGS. 21 and 22, but are repeated in order to allow comparison withFIGS. 27 and 28. FIGS. 27 and 28 show a detail of the plots of FIGS. 25and 26 between 1 ps and 6 ps.

In FIG. 27, the diseased tissue trace lies well below that of the tracesfor healthy tissue and inflamed tissue which follow a similar path.However, in FIG. 28, where E(t)/E_(min) is plotted, a clear distinctioncan be seen between the three traces, the diseases tissue trace beingthe lowest trace, the healthy tissue trace being the middle trace andthe upper trace corresponding to inflamed tissue.

Thus, an image formed by plotting E(t=3.7 ps) can be used to distinguishbetween diseased tissue on the one hand and healthy and inflamed tissueon the other, whereas an image formed by plotting E(t)/E(min) wheret=3.7 ps can be used to distinguish between diseased tissue, healthytissue and inflamed tissue.

When it is desirable to distinguish between three different tissuetypes, it is necessary to chose a time value which emphasises thedifference between all three tissue types and not just two tissue types.

FIG. 29 illustrates two spectra. The upper spectra is produced bysubtracting the spectra formed by calculating E(t)/E_(min) for healthytissue from the spectra formed by calculating E(t)/E_(min) for diseasedtissue. The lower spectra is produced by subtracting the spectra formedby calculating E(t)/E_(min) for healthy tissue from the spectra formedby calculating E(t)/E_(min) for inflamed tissue.

In order to produce an image which clearly shows contrast betweenhealthy and diseased tissue, a time value should be chosen where themagnitude of the upper spectra is greatest. To produce an image whichclearly shows contrast between healthy and inflamed tissue, a time valueshould be chosen where the magnitude of the lower spectra is greatest.However, in order to distinguish between diseased and inflamed tissue,it is necessary to choose points where the difference between the twospectra is also large.

For example, in FIG. 29, both traces show a large negative value atabout 0.7 ps. However, although taking this value of t will produce animage where there is good contrast between healthy tissue and diseasedtissue and between healthy tissue and inflamed tissue, such an image isunlikely to show contrast between diseased tissue and inflamed tissue.However, producing an image of E(t)/E_(min) for t=3.7 ps or t=1.2 pswill show good contrast between all three tissue types.

FIGS. 30 and 31 illustrate images produced by plotting E(t)/E_(min) foreach pixel for t=1.2 ps (FIG. 30) and t=3.7 ps (FIG. 31).

1. A method of imaging a sample, the method comprising: (a) irradiatinga sample with a pulse of electromagnetic radiation, said pulse having aplurality of frequencies in the range from 25 GHz to 100 THz; (b)determining a first parameter at least related to the amplitude of theradiation, which is either reflected from and/or transmitted by thesample, in the time domain; (c) calculating the value of the firstparameter at a first time value relative to the value of the firstparameter at a second time value, wherein the second time valuecoincides with a physical feature of the first parameter with respect totime; and (d) generating an image of the sample by plotting the value ofthe first parameter calculated in step (c) for different spatial pointsover an area of the sample.
 2. A method according to claim 1, whereinsaid physical feature is a predetermined maxima or minima of the firstparameter with respect to time.
 3. A method according to claim 1,wherein step (c) comprises the step of subtracting the value of thefirst parameter at the second time value from the value of the firstparameter at the first time value and step (d) comprises plotting thevalue of the first parameter at the second time value subtracted fromthe value of the first parameter at the first time value for differentspatial points over an area of the sample.
 4. A method according toclaim 1, wherein in step (c), the first time value is one of the values,ascending or descending, to or from, a local minima or maxima and saidsecond time value coincides with said local minima or maxima, andwherein the first time value does not coincide with a local maxima orminima.
 5. A method according to claim 4, wherein step (c) comprises thestep of choosing the second value to correspond to a predeterminedminima of the first parameter and the first time value is chosen fromthe subsequent time values ascending from the said minima.
 6. A methodaccording to claim 1, wherein the first time value is fixed in time foreach point of the sample used to create a single image.
 7. A methodaccording to claim 6, further comprising the step of determining thefirst parameter for a plurality of first time values relative to thefirst parameter for the second time value for each point, the methodfurther comprising the step of generating a plurality of images of thesample corresponding to each of the said first time values.
 8. A methodaccording to claim 7, further comprising the step of selecting theoptimum first time value by selecting the image which provides the bestcontrast.
 9. A method according to claim 1, wherein the first time valueis a fixed time interval from the second time value for each part of thesample used to create a single image.
 10. A method according to claim 9,further comprising the step of determining the first parameter for aplurality of different time intervals relative to the second time valuefor each part, the method further comprising the step of generating aplurality of images of the sample corresponding to each of differenttime intervals.
 11. A method according to claim 10, further comprisingthe step of selecting the optimum time interval by selecting the imagewhich provides the best contrast.
 12. A method according to claim 1,further comprising the steps of: (i) selecting n regions of said sample,where n is an integer of at least 2; (ii) producing n time domainspectra of a second parameter which is at least related to the amplitudeof the radiation reflected from and/or transmitted by the sample, wherethe nth spectra is obtained by averaging the second parameter across thenth region and plotting this averaged value, as it varies with delaytime; (iii) determining the first time value by comparing at least twoof the n spectra derived in step (ii).
 13. A method according to claim1, wherein the first parameter is the amplitude of the radiation.
 14. Amethod according to claim 1, wherein step (a) comprises: irradiating thesample through a member which is transparent to the irradiatingradiation and which abuts the sample, the method further comprising thestep of: obtaining a baseline signal by irradiating the member in theabsence of the sample and detecting the amplitude of the reflectedradiation, and wherein in step (b): the first parameter is determined bysubtracting the baseline signal from the detected amplitude of theradiation reflected from both the member and the sample.
 15. A methodaccording to claim 14, wherein the baseline signal is subtracted fromthe detected reflected radiation in the time domain.
 16. A methodaccording to claim 14, wherein step (b) comprises determining the firstparameter by subtracting the baseline signal from the detected radiationfrom the sample and multiplying the baseline subtracted sample signal inthe frequency domain by the complex Fourier transform of the functionF(t), where F(t) is a non-zero function whose integral between timelimits t_(a) and t_(b) is zero, and where t_(a) and t_(b) are chosen toencompass the part of interest of the reflected or transmitted pulse.17. A method according to claim 1, wherein step (a) comprises:irradiating the sample through a member which is transparent to theirradiating radiation and which abuts the sample, the method furthercomprises: obtaining a baseline signal by irradiating the member in theabsence of the sample and detecting the amplitude of the reflectedradiation, and replacing the sample with a reference object of knownreflectance and measuring a reference signal by measuring the amplitudeof radiation reflected from the reference object, wherein in step (b)comprises: determining the first parameter by subtracting the baselinesignal from the detected radiation and dividing the result by thereference signal in the frequency domain.
 18. A method according toclaim 17, wherein step (b) comprises determining the first parameter bysubtracting the baseline signal from the amplitude of the detectedradiation and dividing the result by the reference signal in thefrequency domain and multiplying by the complex Fourier transform of thefunction F(t) in the frequency domain, where F(t) is a non-zero functionwhose integral between time limits t_(a) and t_(b) is zero, and wheret_(a) and t_(b) are chosen to encompass the part of interest of thereflected or transmitted pulse.
 19. A method according to claim 1,further comprising the step of: replacing the sample with a referenceobject of known reflectance and measuring a reference signal bymeasuring the amplitude of radiation reflected from the referenceobject, and wherein step (b) comprises: determining the first parameterby dividing the detected amplitude with the reference signal in thefrequency domain.
 20. A method according to claim 19, wherein step (b)comprises determining the first parameter by dividing the amplitude ofthe detected radiation from the sample by the reference signal in thefrequency domain and multiplying the result by the complex Fouriertransform of F(t) in the frequency domain, where F(t) is a non-zerofunction whose integral between time limits t_(a) and t_(b) is zero, andwhere t_(a) and t_(b) are chosen to encompass the part of interest ofthe reflected or transmitted pulse.
 21. A method according to claim 19,where:${F(t)} = {❘{\frac{2}{\pi}\left\{ {\frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\alpha})}^{2}}}{\alpha} - \frac{1}{T}} \right\}}}$and wherein α is a constant which is substantially equal to the shortestpulse length of the beam and T is substantially equal to the time whichit takes the beam of radiation to penetrate to the deepest point ofinterest in the sample.
 22. A method according to claim 1, wherein step(b) further comprises: determining the first parameter by multiplyingthe amplitude of the detected radiation in the frequency domain by thecomplex Fourier Transform of function F(t), where F(t) is a non-zerofunction whose integral between time limits t_(a) and t_(b) is zero, andwhere t_(a) and t_(b) are chosen to encompass the part of interest ofthe reflected or transmitted radiation pulse.
 23. A method according toclaim 22, where${F(t)} = {❘{\frac{2}{\pi}\left\{ {\frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\alpha})}^{2}}}{\alpha} - \frac{{\mathbb{e}}^{{- 2}{(\frac{t}{\beta})}^{2}}}{\beta}} \right\}}}$and α and γ are constants.
 24. A method according to claim 22, wherein ais substantially equal to the shortest pulse length of the beam ofpulsed radiation and b is set to be longer than the pulse length.
 25. Anapparatus according to claim 24, further comprising: (i) means forselecting n regions of said sample, where n is an integer of at least 2;(ii) means for producing n time domain spectra of a second parameterwhich is at least related to the amplitude of the radiation reflectedfrom and/or transmitted by the sample, where the nth spectra is obtainedby averaging the second parameter across the nth region and plottingthis averaged value as it varies with delay time; and (iii) meansdetermining the first time value by comparing at least two of the nspectra.
 26. A method of imaging a sample, the method comprising: (a)irradiating a sample with a pulse of electromagnetic radiation, saidpulse having a plurality of frequencies in the range from 25 GHz to 100THz; (b) determining a first parameter at least related to the amplitudeof the radiation, which is either reflected from and/or transmitted bythe sample, in the time domain; (c) calculating the value of the firstparameter at a first time value relative to the value of the firstparameter at a second time value, wherein the second time valuecoincides with a physical feature of the first parameter with respect totime, by dividing the value of the first parameter at the first timevalue with the value of the first parameter at the second time value;and (d) generating an image by plotting the value calculated in step (c)for different points of the sample.
 27. A method of imaging a sample,the method comprising: (a) irradiating a sample with a pulse ofelectromagnetic radiation, said pulse having a plurality of frequenciesin the range from 25 GHz to 100 THz; (b) determining a first parameterat least related to the amplitude of the radiation, which is eitherreflected from and/or transmitted by the sample, in the time domain; (c)calculating the value of the first parameter at a first time valuerelative to the value of the first parameter at a second time value; and(d) generating an image of the sample by plotting the value of the firstparameter calculated in step (c) for different spatial points over anarea of the sample, wherein the time interval between the first andsecond time values is constant for all points used to generate theimage.
 28. method of imaging a sample, the method comprising: (a)irradiating a sample with a pulse of electromagnetic radiation, saidpulse having a plurality of frequencies in the range from 25 GHz to 100THz; (b) determining a first parameter at least related to the amplitudeof the radiation, which is either reflected from and/or transmitted bythe sample, in the time domain; and (c) generating an image of thesample by plotting the value of the first parameter at a first timevalue, said first time value being determined by: (i) selecting nregions of said sample, where n is an integer of at least 2, (ii)producing n time domain spectra of a second parameter which is at leastrelated to the amplitude of the radiation reflected from and/ortransmitted by the sample, where the nth spectra is obtained byaveraging the second parameter across the nth region and plotting thisaveraged value, as it varies with delay time; (iii) determining thefirst time value by comparing at least two of the n spectra derived instep (ii).
 29. A method according to claim 28, wherein step (ii)comprises generating a preliminary image of an area of the sample byplotting a preliminary imaging parameter which is at least related tothe amplitude of the radiation reflected from and/or transmitted by thesample; and selecting said n regions from said preliminary image.
 30. Amethod according to claim 29, wherein said n regions are selected suchthat said n regions being selected such that the level of thepreliminary imaging parameter within each region is substantiallyconstant and the level of the preliminary imaging parameter is differentbetween at least two of the said regions.
 31. A method according toeither of claim 29, wherein the preliminary imaging parameter is thevalue of the first and/or second parameter at a time which correspondsto a physical feature of the time domain spectra of the first and/orsecond parameter.
 32. A method according to claim 29, where n=3 and thevalue of the preliminary imaging parameter differs between said threeselected regions, the first time value being chosen to coincide with aregion where there is a difference between the three spectra.
 33. Amethod according to claim 28, wherein each spectra is normalised to thevalue of a third parameter at a time value which coincides with aphysical feature of the second parameter with respect to time, prior tocomparing the spectra.
 34. A method according to claim 28, wherein step(iii) comprises performing a mathematical operation to combine at leasttwo of the spectra.
 35. A method according to claim 34, wherein thedifference between at least two of the spectra is calculated.
 36. Amethod according to claim 34, wherein at least one of the spectra isdivided by another one of said n spectra.
 37. A method according toclaim 28, wherein n=2.
 38. A method according to claim 28, wherein thesecond parameter and the first parameter are the same.
 39. An apparatusfor imaging a sample, the apparatus comprising: a source for irradiatinga sample with a pulse of electro-magnetic radiation, said pulse having aplurality of frequencies in the range from 25 GHz to 100 Thz; detectormeans for detecting the amplitude of radiation reflected from and/ortransmitted by the sample; means for determining a first parameter atleast related to the amplitude of the radiation in the time domain;calculating means configured to calculate the value of the firstparameter at a first time value with respect to the value of the firstparameter at a second time value, wherein the second time valuecoincides with a physical feature of the dataset of the first parameterwith respect to time, and imaging means for generating an image of thesample by plotting the value of the first parameter calculated by thecalculating means for different spatial points over an area of thesample.
 40. An apparatus according to claim 39, wherein the calculatingmeans is configured to calculate the value of the first parameter forthe first time value with respect to the second time value for aplurality of different first time values for each point of the sampleused to generate the image.
 41. An apparatus according to claim 40,further comprising means to produce a plurality of images for differentfirst time values.
 42. An apparatus according to claim 41, wherein theapparatus comprises means to scan through the plurality of images. 43.An apparatus according to claim 39, wherein the first time value isconstant for all points used to generate a single image.
 44. Anapparatus according to claim 43, further comprising means to optimisecontrast in the image by selecting the first time value which providesthe best contrast.
 45. An apparatus according to claim 39, wherein thetime interval between the first time value and the second time value isconstant for all points used to generate a single image.
 46. Anapparatus according to claim 45, further comprising means to optimisecontrast in an image by selecting the fixed time interval which providesthe best contrast.
 47. An apparatus for imaging a sample, the apparatuscomprising: a source for irradiating a sample with a pulse ofelectro-magnetic radiation, said pulse having a plurality of frequenciesin the range from 25 GHz to 100 THz; detector means for detecting theamplitude of radiation reflected from and/or transmitted by the sample;means for determining a first parameter at least related to theamplitude of the radiation in the time domain; calculating meansconfigured to calculate the value of the first parameter at a first timevalue with respect to the value of the first parameter at a second timevalue, and imaging means for generating an image of the sample byplotting the value of the first parameter calculated by the calculatingmeans for different spatial points over an area of the sample, whereinthe time interval between the first and second time intervals isconstant for all points used to generate the image.
 48. An apparatus forimaging a sample, the apparatus comprising: (a) means for irradiating asample with a pulse of electromagnetic radiation, said pulse having aplurality of frequencies in the range from 25 GHz to 100 THz; (b) meansfor determining a first parameter, at least related to the amplitude ofthe radiation, which is either reflected from and/or transmitted by thesample, in the time domain; and (c) means for generating an image of thesample by using the value of the first parameter at a first time value,said apparatus being configured to determine said first time value by:(i) selecting n regions of said sample, where n is an integer of atleast 2; (ii) producing n spectra in the time domain, where the nthspectra is obtained by plotting a second parameter which is at leastrelated to the amplitude of the radiation reflected from and/ortransmitted by the sample, averaged across the nth area, against delaytime; iii) determining the first time value by comparing at least two ofthe n spectra.